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Abstract 

Philaenus spumarius is a widespread insect species in the Holarctic region. Here, by focusing on the mtDNA gene COI but 
also using the COM and Cyt b genes and the nuclear gene EF-1a, we tried to explain how and when its current 
biogeographic pattern evolved by providing time estimates of the main demographic and evolutionary events and 
investigating its colonization patterns in and out of Eurasia. Evidence of recent divergence and expansion events at less 
than 0.5 Ma ago indicate that climate fluctuations in the Mid-Late Pleistocene were important in shaping the current 
phylogeographic pattern of the species. Data support a first split and differentiation of P. spumarius into two main 
mitochondrial lineages: the "western", in the Mediterranean region and the "eastern", in Anatolia/Caucasus. It also supports 
a following differentiation of the "western" lineage into two sub-lineages: the "western-Mediterranean", in Iberia and the 
"eastern-Mediterranean" in the Balkans. The recent pattern seems to result from postglacial range expansion from Iberia 
and Caucasus/Anatolia, thus not following one of the four common paradigms. Unexpected patterns of recent gene-flow 
events between Mediterranean peninsulas, a close relationship between Iberia and North Africa, as well as high levels of 
genetic diversity being maintained in northern Europe were found. The mitochondrial pattern does not exactly match to the 
nuclear pattern suggesting that the current biogeographic pattern of P. spumarius may be the result of both secondary 
admixture and incomplete lineage sorting. The hypothesis of recent colonization of North America from both western and 
northern Europe is corroborated by our data and probably resulted from accidental human translocations. A probable 
British origin for the populations of the Azores and New Zealand was revealed, however, for the Azores the distribution of 
populations in high altitude native forests is somewhat puzzling and may imply a natural colonization of the archipelago. 
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Introduction 

Distribution patterns of animals and plants have faced dramatic 
changes throughout time and are influenced by ecological 
requirements and historical factors. In the northern hemisphere, 
Quaternary long-term glacial (cold) and interglacial (warm) 
climatic cycles that started about 2.6 million years (Ma) ago [1] 
have strongly influenced the species distributions and range sizes 
and, as a consequence, have affected the genetic structure of their 
populations [2,3]. Evidence from numerous studies suggests that 
southern European regions of Iberia, Italy and the Balkans and 



areas near the Caucasus and western Asia, acted as glacial refugia 
for temperate species during cold periods [4,5,6]. Recent work 
indicates that temperate refugia might not have been restricted to 
the three southern peninsulas and that cryptic northern refugia 
might have existed in central, western, eastern and even northern 
Europe in the Late Pleistocene [7,8,9]. The relative impact of the 
post-glacial colonization history and more recent processes such as 
gene flow and population fluctuations, strongly depend on the 
dispersal mode and ability of the species [10,11]. 

Genetic analyses have proven to be useful for a more detailed 
understanding of post-glacial expansions of several animals and 
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Figure 1 . Sampling locations of Philaenus spumarius in (a) Europe and Anatolia (b) North America and (c) New Zealand in each 
geographic region. 1 - Azores; 2 - Iberian Peninsula; 3 - Morocco; 4 - France; 5 - United Kingdom; 6 - Belgium; 7 - Italian Peninsula; 8 - Sicily; 9 - 
Slovenia; 1 0 - Balkans (Bulgaria; Greece; European Turkey); 1 1 - Anatolian Peninsula; 1 2 - Finland; 1 3 - North America (Canada and United States of 
America); 14 - New Zealand. Circle sizes are proportional to the number of individuals. Circles: green - "western-Mediterranean" mtDNA group; red - 
"eastern-Mediterranean" mtDNA group; blue - "eastern" mtDNA group. Circle sizes are proportional to the number of samples. 
doi:10.1371/journal.pone.0098375.g001 



plants [3,12]. Mitochondrial DNA (mtDNA), due to its particular 
characteristics, has been widely used in determining population 
dynamics and phylogeographic divergence in recent times, such as 
the Quaternary period [13]. Nevertheless, the signal of deeper 
history can be obscured by homoplasy or saturation resulting from 
high mutation rate. On the other hand, reconstructing evolution- 
ary histories using individual genes (gene trees) could lead to 
misrepresentation of population or species histories because in this 
case mtDNA, which reflects matrilineal history, might not 
represent the overall lineage history of the species. Also, if multiple 
population divergences or speciation events were closely spaced in 
time, a single gene tree might be 'incorrect' by chance due to the 
random nature of lineage sorting during the coalescence process 
[14]. Therefore, the use of multiple types of molecular markers is 
recommended. 

Insects have been widely used as models for animal biogeo- 
graphical studies (e.g., [15,16,17]). The meadow spittlebug 
Philaenus spumarius (Linnaeus, 1758) (Hemiptera, Aphrophoridae) 
is a widely investigated species, very suitable for genetic and 
ecological studies. It is a highly polyphagous insect which can be 
found in a variety of terrestrial plant communities and habitats, 
being the most common species within the genus Philaenus [18,19]. 
It is widespread across the Palaearctic region from where it is 
native [20] having also colonized the Azores [21,22] and has been 
introduced in the Nearctic region [20] and New Zealand [23]. 
The meadow spittlebug is very sensitive to humidity and 
temperature, especially in the earlier stages of its life cycle, which 
limits its range [24] . A remarkable example was reported for some 
North American populations where a northward range shift, 
probably as a result of climatic changes, was detected by [25]. This 
species shows a well studied dorsal colour polymorphism with 
eleven main described phenotypes which can be divided in 
melanic and non-melanic forms [20] . The phenotype frequencies 
differ among populations, probably due to the effects of natural 
selection under different habitats, climatic conditions and preda- 
tion pressure (reviewed in [18,19]). Recent studies on the genetic 



diversity of P. spumarius have given insights on its evolutionary 
history suggesting two routes of post-glacial colonization of higher 
latitudes in Europe and indicating a probable western European 
origin for North American populations [26,27]. 

In the present study we tried to explain how and when the 
current biogeographic pattern of P. spumarius evolved by (i) 
providing time estimates of the main demographic and evolution- 
ary events with focus on the populations occurring in the main 
Mediterranean peninsulas; and, (ii) investigating the colonization 
patterns out of Eurasia, namely of north-western Africa, North 
America, and the islands of the Azores and New Zealand. 

Material and Methods 

Ethics Statement 

The field sampling was carried out on private lands with 
owners' permissions. The studied species, Philaenus spumarius, is 
considered a widespread species across the Palaearctic and the 
Nearctic regions, being a crop pest in some locations of USA and 
Canada. It is not an endangered or protected species. 

Sampling 

A total of 196 specimens of P. spumarius were collected or sent by 
collaborators between 2007 and 201 1 from 75 sampling locations 
across Europe, two from Anatolia, five from North Africa, three 
from North America and one from New Zealand (Fig. 1 and Table 
SI). Adult insects were captured using a sweep net suitable for low- 
growing vegetation and an entomological aspirator. In some cases, 
larval stages were collected by hand. Specimens were preserved in 
absolute ethanol or dried in silica gel and stored at room 
temperature. 

DNA extraction, amplification and sequencing 

Entire larval stage specimens were used for DNA extraction 
while in the adults the wings and abdomen were removed and only 
the thorax and head were used. Genomic DNA was extracted 
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using the E.Z.N.A. Tissue DNA Isolation kit (Omega Bio-Tek) and 
a 800 bp fragment of the 3'-end of the mitochondrial gene 
cytochrome c oxidase subunit I (COI) was amplified by polymerase 
chain reaction (PCR) using the primers Cl-J-2195 (5'— 
TTGATTTTTTGGTCATCCAGAAGT-3') and TL2-N-3014 
(5 '-TCCAATGCACTAATCTGCCATATTA-3') [28]. Due to 
DNA degradation in the New Zealand samples, a new set of 
primers, COI-PspF (5 '-GTATAGATGTTGATAC AC GTGC- 
3') and COI-PspR (5'-TC C AGTAAATAAAGGGTATC-3 ') 
was designed to amplify an informative smaller fragment with 
300 bp of COI that included the variable sites that differentiate the 
different haplogroups. Fragments with 500 bp of the mitochon- 
drial genes cytochrome c oxidase subunit II (COII) and 
cytochrome b (cyt b) were also amplified using the primers TL2- 
J3033 (5'-GATATGGCAGAAATAGTGCA-3') and C2-N3665 
(5'-CCACAAATTTCTGAACACTG-3') and CB-N3665 (5'- 
GTCCTACCATGAGGTCAAATATC-3') and CB-N11526 
(5'-TTCAACTGGTCGTGCTCC-3'), respectively [29]. The 
three mitochondrial genes were sequenced initially in a subset of 
samples and similar genetic patterns and level of polymorphism 
were observed for the three mitochondrial genes on a preliminary 
analysis (results not shown). Thus, cytochrome c oxidase subunit I 
gene was sequenced for all individuals and cyt b and COII were 
only sequenced in a representative subset. The nuclear gene 
elongation factor- la (EF-la) is widely used in insect genetic 
studies. Therefore, we chose to sequence it in a subset of 
individuals that covered all the geographical areas of the study. A 
700 bp fragment of the nuclear gene EF-la was amplified using 
the primers M3 (5 '-C AC ATYAAC ATTGTC GTS ATYGG— 3 ') 
and rcM44.9 (5'-CTTGATGAAATCYCTGTGTCC-3') [30]. 
For COI gene, PCR was performed in a 12.5 uL reaction volume 
containing: 1 |jM of each primer, 0. 1 mM dNTPs, 1 mM MgCl 2 , 
2.5 uL 5x Colorless GoTaq Flexi Buffer, 0.02U GoTaq DNA 
Polimerase (Promega) and approximately 30 ng of DNA. The 
PCR conditions were: an initial denaturation step at 94 °C for 
3 min, followed by 35 cycles of denaturation at 94 °C for 30 sec, 
annealing at 50 °C for 45 sec and extension at 72 °C for 1 min, 
with a final extension period at 72 °C for 7 min. The same PCR 
conditions were used for COII and cyt b genes except for 
annealing temperature where a touch up between 52.5 °C and 56 
°C for COII and between 47 °C and 54 °C for cyt b was 
performed. Nuclear EF-la gene PCR was performed in a 20 uL 
reaction volume containing: 0.6 uM of each primer, 0.2 mM 
dNTPs, 1.125 mM MgCl 2 , 0.8 uL BSA (10 ug/mL), 4.0 uL 5x 
Colorless GoTaq Flexi Buffer, 0.05U GoTaq DNA Polimerase 
(Promega) and approximately 30 ng of DNA. PCR conditions 
used were: an initial denaturation step at 95 °C for 5 min, followed 
by 40 cycles of denaturation at 95 °C for 45 sec, annealing at 59 
°C for 35 sec and extension at 72 °C for 1 min, with a final 
extension period at 72 °C for 10 min. All PCR products were 
purified with SureClean (Bioline) following the manufacturer's 
protocol, sequenced using the forward and the reverse primers 
with the BigDye Terminator v3. 1 Cycle Sequencing Kit (Applied 
Biosystems) and analysed on a genetic analyser ABI PRISM 310 
(Applied Biosystems). 

Molecular and population structure analyses 

Sequences were verified and edited using the software 
Sequencher v. 4.0.5 (Gene Codes Corporation) and BioEdit v. 
7.0.9 [31]. They were then aligned using Mafet v. 7.029b (http:// 
mafft.cbrc.jp/alignment/software/) and converted in the appro- 
priate format with Concatenator v. 1.1.0 [32]. For nuclear EF- 
la sequences, haplotype phase from heterozygous individuals for 
base positions and length-variable regions was determined using 



Champuru v. 1.0 [33]. Phylogenetic analysis using the Maximum 
Parsimony (MP), Maximum Likelihood (ML) and the Bayesian 
inference (BI) methods were performed for concatenated mtDNA 
genes (COI, COII and cyt b) and for the nuclear gene EF-la in 
Paup v. 4.0.d99 [34] and in MrBayes v. 3.1.2 [35]. For MP and 
ML analysis a heuristic search was performed using 1 00 replicates 
and branch support was obtained by performing 1000 replicates of 
non-parametric bootstrap. Gaps were treated as a fifth base in MP. 
The BI analysis was performed using the Monte Carlo Markov 
Chain (MCMC) method iterated for 2 000 000 generations, with a 
sampling frequency of 1500 generations and a burn-in of 1000. 
For each dataset the best fit model of sequence evolution was 
estimated using Modeltest v. 3.7 [36] under the Akaike 
information criterion (AIC). Elongation factor- la sequences of P. 
spumarius and P. italosignus from [26] and available at NCBI 
Genbank were added to our nuclear matrix and included in the 
phylogenetic analysis (GenBank accession numbers: JF309079 and 
JF309081-JF309095). Philaenus italosignus was used as outgroup in 
all phylogenetic analysis. Polymorphic sites and mtDNA haplo- 
types for COI, COII and cyt b genes were calculated using Mega 
v. 5.0 [37] and a median-joining haplotype network was 
constructed using Network v. 4.5.0.1 (Fluxus Technology Ltd. 
2004). For COI mtDNA gene, haplotype (h) and nucleotide 
diversities (n) were calculated for each geographical region 
(defined as numbers in Fig. 1) and an analysis of molecular 
variance (AMOVA) was performed using Arlequin v. 3.5 [38] to 
assess population genetic structure of P. spumarius. The groupings 
were based in the several sub-regions of Europe, America, Africa 
and Asia. This analysis produces estimates of variance components 
and F-statistic analogues, designated as CD-statistics, reflecting the 
correlation of haplotypes at different levels of hierarchical 
subdivision. Groupings with the highest significant <1>ct value in 
AMOVA should reflect the most probable geographical subdivi- 
sions [39]. 

Divergence time estimates 

We used the software package Beast v. 1.7.0 [40] and the 
mtDNA gene COI to estimate divergence times of nodes of 
interest, as well as their demographic history via Bayesian Skyline 
plots (BSPs). For each dataset the best fit model of sequence 
evolution was estimated using Modeltest v. 3.7 under the Akaike 
information criterion (AIC) and a piece-wise constant Bayesian 
skyline tree prior was selected with 10 groups. Two additional 
analyses of our data with 5 and 15 groups to assess the impact of 
the number of groups on the final result were conducted. These 
analyses did not reveal a significant impact on the overall result, 
whether on the shape of the Bayesian Skyline plot or on the 
estimation of divergence times. Preliminary runs using the 
uncorrelated lognormal relaxed clock revealed a posterior 
distribution of the c? r ("CoefficientOfVariation") parameter 
consistently abutting 0, suggesting that the COI partitions do 
not significantly deviate from a strict clock assumption. Therefore, 
we employed a strict molecular clock for each dataset with a 
normal prior distribution on the substitution rate with a mean of 
0.0354, following the conserved rate of 3.54% per million years as 
suggested by [41], and a standard deviation of 0.005 to account for 
rate uncertainty. Two Markov chains Monte Carlo (MCMC) of 50 
000 000 generations, sampled at every 5000 th iteration, were 
conducted and combined with LogCombiner v. 1.6.1 [40]. 
Tracer v. 1.4 [42] was used to assess the convergence and 
mixing for all model parameters and to create the Bayesian Skyline 
plots. 
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Demographic analyses and neutrality tests 

Neutrality tests of Tajima's D [43] and Fu's F statistics [44] were 
performed using Arlequin. These statistics are widely used with 
molecular data to detect changes in population size and/or 
estimating deviations from neutrality, assuming a constant 
population size at mutation-drift equilibrium. Thus, significant 
negative values of Tajima's D and Fu's Fs are considered to be 
evidence of expanding populations. Signatures of population 
expansion can also be detected through the frequency distribution 
of the number of pairwise differences between haplotypes and thus 
statistics based on the mismatch distribution and taking into 
account the Sudden Expansion Model [45] were also performed to 
detect and estimate the time of population growth. Estimated 
expansion values were obtained using Arlequin and graphics of 
frequency distribution using DnaSP v. 5 [46] . To test the observed 
mismatch distribution goodness-of-fit to the Sudden Expansion 
Model and to obtain confidence intervals around the estimated 
mode of mismatch distribution, 1 000 permutation replicates were 
used [47]. Statistically significant differences between observed 
and expected distributions were evaluated with the sum of the 
square deviations (SSD) and Harpending's raggedness index (fig) 
[48,49]. 

Timing of the demographic expansion as well as the 95% 
confidence interval for each mitochondrial haplogroup was 
estimated by converting the expansion time parameter T, 
generated by Arlequin, to time (t) in years using the formula x 
= 2ut, where u is the mutation rate per nucleotide per year 
multiplied by sequence length (i.e. number of nucleotides), and t is 
the time since population expansion in years [45,49]. We assumed 
a generation time of one year [24] and the conserved evolutionary 
rate of 3.54% per million years suggested by [41] for insect 
mitochondrial gene COI. 

Results 

A total fragment of 539 bp was obtained for the mitochondrial 
gene COI in 190 samples. The remaining samples from Finland, 
Turkey, Canada (GenBank accession numbers: KJ699232— 
KJ699234) and New Zealand were not included in the main 
analysis due to their reduced size. In the 190 individuals there 
were a total of 71 haplotypes (GenBank accession numbers: 
KC 1 1 1 886 - KC 1 1 1 956) of which 44 occurred only once (Table 
S2). Of a total of 539 sites sequenced, 53 were polymorphic but 
only 26 were parsimony informative. For the mtDNA COII 
(495bp) and cyt b (434bp) genes, 47 individuals were sequenced 
and a total of 14 (GenBank accession numbers: KF280589 - 
KF280602) and 18 haplotypes (GenBank accession numbers: 
KF280603 - KF280620) were found, respectively. As commonly 
observed for insects [28], nucleotide sequences were A+T rich 
(approximately 7 1 %). No gaps or early stop codons were detected 
in the 3 mtDNA genes sequences suggesting that all of them are 
functional mitochondrial DNA copies. 

From the 24 individuals sequenced, we were able to successfully 
sequence a fragment of the nuclear gene EF-la for only 13 
individuals of P. spumarius. Almost all sequences exhibited double 
peaks due to the frameshift resulting from indels (length-variable 
regions) located at several sites of the intron in the EF-la 
gene. Ten individuals were heterozygous in respect to indels and/ 
or to base positions and the phased haplotypes (alleles) were 
differentiated by adding the letter a or b at the end of 
name (Table SI) (GenBank accession numbers: KF280621 
- KF280642). 



Phylogenetic and population structure analyses 

Phylogenetic trees obtained for concatenated mtDNA genes and 
for a subset of P. spumarius individuals by the three methods, MP, 
ML and BI, presented a congruent topology. Maximum likelihood 
(Fig. 2), MP (not shown) and BI (not shown) trees, revealed the 
existence of two main haplotype groups: the "western" and the 
"eastern". The "western" is divided in the "western-Mediterra- 
nean" and the "eastern-Mediterranean" sub-groups. The same 
phylogeographic pattern was found in the COI median-joining 
haplotype network (Fig. 3) and also in the COII and Cyt b median- 
joining haplotype networks (Fig. SI and S2). The "eastern" 
haplogroup includes haplotypes from a wide geographical area, 
including northern Anatolia (Cerkes), Finland, Belgium, the UK 
(Aberdare - Wales) and eastern North America (New Hampshire - 
USA). In the "western-Mediterranean" group, the most common 
haplotype (H29) and several derived haplotypes, differing by one 
or two mutational steps, are shared between populations from the 
Iberian Peninsula, Morocco, France, Belgium, Italian Peninsula, 
Sicily and one individual from Balkans (HI 8). A group of 
haplotypes derived from H29 (H23, H24, H25, H70 and H71) 
includes samples from the Azores, western North America (British 
Columbia - Canada), eastern North America (Michigan - USA) 
and the UK, differing by two or three mutations. In the "eastern- 
Mediterranean" group, a similar star-like pattern is present with 
rare haplotypes connected to the most common (H57), usually by 
one mutational step. This group encompasses populations from 
the Balkans (Greece, Bulgaria and European Turkey), Slovenia, 
Italian Peninsula and Sicily. This lineage is also present in five 
samples from the eastern part of the Iberian Peninsula (H56) 
(Figs. 1 and 2). 

The haplotype median-joining network based on smaller 
sequences of 289 bp of the COI gene, used in this analysis to 
include the samples from New Zealand, showed a total of 31 
haplotypes and a pattern congruent to that observed for the 539 
bp of COI region, with the same three distinct groups ("western- 
Mediterranean", "eastern-Mediterranean" and "eastern"). The 
haplotype H14 belonging to the "western" haplogroup was found 
in the three New Zealand individuals and this same haplotype was 
shared with seven individuals from the UK, five from the Azores, 
three from Canada and four from the USA (Fig. S3). 

In the MP (Fig. 4) and BI (not shown) phylogenetic trees 
obtained for nuclear gene EF-la three main groups could also be 
distinguished: elade A, clade B and clade C. However, these groups 
were not totally congruent to the groups found for the 
concatenated mtDNA genes. The clade A includes samples from 
our "eastern-Mediterranean" haplogroup and also from Georgia, 
Bulgaria, Hungary, Greece and Italy (E2 clade, in [26]). The clade B 
includes individuals from our "western-Mediterranean" group and 
from Portugal, Spain and Italy (E3 clade, in [26]). The clade C, 
however, is constituted by individuals from our three main 
mitochondrial groups and from Russia, Norway, Alps, Crimea, 
Poland and Ukraine (El clade, in [26]). Although the three groups 
have good bootstrap support, the branching order is unsolved 
since there is very low support for the branch clustering the clades A 
and B (Fig. 4). We also observed that both alleles of the most 
heterozygous samples are clustered within the same clade (Clade C - 
Arrabida_5, Spain_2, UK_7, Belgium_l, Slovenia_l.l and 
USA_2 and Clade B - Morocco_6.1), with the exception for the 
Azores_l, Italy_2 and Kecan_Nl samples which have one allele in 
the clade C and the other allele in the clades A or B. 

Genetic variability and diversity 

Mitochondrial haplotype diversity of COI was generally high (h 
> 0.6000), except for the Azores (h = 0.4000) and Slovenia (h = 
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Figure 2. Maximum Likelihood tree based on the 3 concatenated mtDNA genes (COI, COM and cyt b) (1 527bp). Values above branches 
correspond to MP and ML bootstrap values (only values > 50% are shown) and values below branches correspond to Bayesian posterior probability. 
PT - Portugal; SP - Spain; UK - United Kingdom; BL - Belgium; FN - Finland; SC - Sicily; TK - European Turkey; AT - Anatolia; USA - United States of 
America; CA - Canada. 
doi:10.1371/journal.pone.0098375.g002 



0.0000) (Table 1), that may be related to the low sample size. On 
the other hand, Anatolia and the Balkans have the highest values 
of haplotype diversity (h = 0.9048 and h = 0.9128, respectively). 
Interestingly, Finland and North America also have high 
haplotype diversity (h = 0.8971 and h = 0.8939). Analysing 
nucleotide diversity, Anatolia, United Kingdom and North 



America have the highest values (n = 0.009895, Tt = 0.009318 
and Tt = 0.008236, respectively). The high values of both 
haplotype and nucleotide diversities detected in Anatolia and 
North America are likely a result of the presence of different 
mitochondrial lineages in these regions. The AMOVA performed 
for the P. spumarius groupings defined revealed that most of the 



PLOS ONE | www.plosone.org 



5 



June 2014 | Volume 9 | Issue 6 | e98375 



Phylogeography of the Spittlebug Philaenus spumarius 




Figure 3. Median-joining haplotype network of Philaenus spumarius sampled geographic regions for mitochondrial gene COI 
(539bp). Size of the circles is in proportion to the number of haplotypes. Branches begin in the centre of the circles and their sizes are in proportion 
to the number of mutations. 
doi:1 0.1 371 /journal.pone.0098375.g003 



genetic variation (47% and 43.87%) is explained by genetic 
differences within populations and not by geographic subdivisions 
(Table S3). 

Divergence times 

The mean ages and 95% highest posterior density (HPD) 
determined for the TMRCA of mtDNA COI haplogroups are 
presented in Table S4. Estimated divergence times of all 
mitochondrial groups are less than 0.5 Ma. The "eastern" 
mtDNA group seems to be the oldest having diverged from the 
"western" mtDNA group at approximately 0.270 Ma ago and 
then begun its diversification around 0.190 Ma ago (0.374-0.056 
Ma ago). The split of the "western" mtDNA group into the 
"eastern-Mediterranean" and the "western-Mediterranean" sub- 
groups was estimated to have occurred at approximately 0.146 Ma 
ago (0.243-0.067 Ma ago), while the TMRCA of both groups was 
quite similar and was estimated to be around 0.080 Ma ago. The 
confidence interval associated with our time estimates, however, is 
broad and the mutation rate of 3.54% per million years, in which 
these calculations are based, was estimated for COI in tenebrionid 
beetles [41]. Therefore, our results should be treated with caution 
and regarded as the best approximations given the current 
methods and calibrations [50,51]. 



Demographic analyses and neutrality tests 

The demographic history of P. spumarius populations was 
analysed separately for the three COI groups. The distribution 
of pairwise nucleotide differences (mismatch distribution) showed 
that the "western-Mediterranean" and the "eastern-Mediterra- 
nean" groups exhibited a smooth and unimodal shape while the 
"eastern" group revealed a slightly bimodal curve (Fig. S4). All 
distributions, except for the "eastern" group distribution, were 
consistent with sudden and spatial population expansions. The 
observed raggedness index was low for all groups and both Pssd 
and -Prag showed that the observed distributions did not differ 
significandy from those expected under a sudden and a spatial 
population expansion model (Table 2). Negative significant 
deviations from neutrality were detected with Tajima's D and 
Fu's .F statistics for the "western-Mediterranean" and the "eastern- 
Mediterranean" mtDNA groups, which corroborate the hypoth- 
esis of past population expansion events. The "eastern" group 
presented a non significant Rvalue with Tajima's D test although 
the Fu's F was significant (Table 3) indicating that it may have 
undergone negligible population growth. Demographic recon- 
structions (BSPs) for all mtDNA groups suggest a trend of 
population growth (Fig. S5), with a more evident demographic 
expansion in the "western-Mediterranean" and the "eastern- 



PLOS ONE | www.plosone.org 



6 



June 2014 | Volume 9 | Issue 6 | e98375 



Phylogeography of the Spittlebug Philaenus spumarius 



62/0.86 



53/0.90 



hapEF2 



54/- 



- Arrabida_5a 
98/1 r^ hapEFIl" 
64 H L hapEF3 
■ Azores_1 b 
USA_2a 
L Spain_2a 
r UK_7b 
UK 7a 



hapEF5 
- hapEF12 
hapEF13 
Arrabida_5b 
Turkey_3 
- Belgium_1b 
ltaly_2a 
Slovenia_1.1b 
Kegan_N1b 



54/0.88 




a 



hapEF4 
Belgium_1a 

hapEF6 

Spain_2b 
- Slovenia_1.1a 
USA_2b 
ltaly_2b 

hapEF14 

r Morocco_6.1a 

99/1 iM_ Morocco_6.1b 
97/0.94r Belgium_4 



Western-Mediterranean 
Sub-lineage 
+ 

Eastern-Mediterranean 
Sub-lineage 
+ 

Eastern 
Lineage 
+ 

E1 Clade 



Sintra_10 

99/0.66 1 Azores_1a 



-hapEFI 



Western-Mediterranean 
Sub-lineage 

+ 

E3 Clade 



64/0.99. 



hapEF15 
hapEF8 
hapEF7 

hapEF10 

Kegan_N1a 



hapEF9 



Eastern-Mediterranean 
Sub-lincage 

+ 

E2 Clade 



P. italosignus 



20.0 



Figure 4. Maximum Parsimony tree based on nuclear gene elongation factor- 1 a. Values above branches correspond to MP bootstrap (only 
values > 50% are shown) and Bayesian posterior probability values. Black: GenBank sequences (see [26]); blue individuals correspond to the eastern 
mtDNA group; red individuals correspond to the eastern-Mediterranean mtDNA group and green individuals to the western-Mediterranean mtDNA 
group. 
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Mediterranean" mtDNA groups, and very slight or absent 
population growth for the "eastern" lineage. 

Assuming a slight population growth for the "eastern" group, 
the timing of demographic expansion was estimated to have 
occurred at approximately 0.121 Ma ago (0.195-0.041 Ma ago), 



while for the "western-Mediterranean" and the "eastern-Medi- 
terranean" haplogroups it was more recent, at 0.055 Ma ago 
(0.085-0.021 Ma ago) and 0.058 Ma ago (0.079-0.027 Ma ago), 
respectively. The spatial expansion for the three mtDNA groups 
was estimated to be slightly more recent than demographic 
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Table 1. Number of individuals, number of haplotypes and genetic diversity indices calculated for geographic regions of Philaenus 
spumarius and for mitochondrial gene Cytochrome c oxidase I (COI). 





Geographic regions 


Number of 


Number of 


Haplotype diversity (A) 


Nucleotide diversity (n) 




individuals 


haplotypes 










Morocco 


7 


3 


U.DDD/ T7 


0 1 598 


n nm 41 4 j-/ 

U.UU 1 1 1 +/ 


n nm 

U.UU 1 jjo 


Iberian Peninsula 


63 


19 


0.7798 +/- 


0.0493 


0.003133 +/- 


0.002062 


Azores 


5 


2 


0.4000 +/- 


0.2373 


0.000742 +/- 


0.000944 


Western Europe 


9 


4 


0.7500 +/- 


0.1121 


0.006597 +/- 


0.004194 


Slovenia 


3 


1 


0.0000 +/- 


0.0000 


0.000000 +/- 


0.000000 


Italy 


17 


8 


0.8162 +/- 


0.0815 


0.004666 +/- 


0.002957 


Balkans 


40 


21 


0.91 28 +/- 


0.0303 


0.003825 +/- 


0.002431 


Anatolian Peninsula 


7 


5 


0.9048 +/- 


0.1033 


0.009895 +/- 


0.006223 


United Kingdom 


10 


3 


0.6889 +/- 


0.1038 


0.009318 +/- 


0.005592 


Finland 


17 


9 


0.8971 +/- 


0.0534 


0.006603 +/- 


0.003950 


North America 


12 


8 


0.8939 +/- 


0.0777 


0.008236 +/- 


0.004917 



Western Europe: Belgium and France; Italy: Italian peninsula and Sicily. 
doi:1 0.1 371/journal.pone.0098375.t001 



expansion, at 0.080 Ma ago (0.152-0.035 Ma ago) for the 
"eastern" group, at 0.056 Ma ago (0.075-0.028 Ma ago) for the 
"eastern-Mediterranean" and at 0.054 Ma ago (0.074-0.022 Ma 
ago) for the "western-Mediterranean" group. 

Discussion 

Biogeographical patterns, divergence time and demographic 
events in Philaenus spumarius 

Our time estimates indicate that the evolutionary history of P. 
spumarius is most likely related to climate changes of the Pleistocene 
epoch ( — 2.588—0.01 17 Ma ago [1]). Divergence within species is 
estimated to be recent (no more than 0.5 Ma) occurring most 
probably in the Middle/Late Pleistocene. The biogeographical 
pattern of P. spumarius obtained from mtDNA genes shows the 
differentiation of two main mtDNA lineages, the "western" in the 
Mediterranean region and the "eastern" in Anatolia/Caucasus. 
Within the "western" lineage we observed two sub-lineages: the 
"western-Mediterranean" centred in the Iberian Peninsula and 
the "eastern-Mediterranean" centred in the Balkans (Fig 1). This 
pattern was first found by [27], who gave more emphasis to the 
westernmost lineages, and was later corroborated by [26] that 
brought some insights regarding possible refugia located in eastern 
Europe and western Asia. According to our data, during the 
Mindel or Riss glacial period, a first split between the western and 
eastern populations seems to have occurred with the diversification 
of the "western" lineage in the Mediterranean region, and of the 
"eastern" lineage maybe in Anatolia and surrounding area of the 
Caucasus, or even in territories of western Asia, as also suggested 
by [26]. During the following interglacial, the "eastern" lineage 
seems to have suffered a negligible population growth compared 
with a more significant demographic expansion of the "western" 
lineage, which appears to have later retracted to two Mediterra- 
nean refugia, the Iberian Peninsula and the Balkans, where it 
diverged into two sub-lineages (the "western-Mediterranean" and 
the "eastern-Mediterranean"), maybe during the Wurm glacial. 
After that period, the "eastern-Mediterranean" lineage centred in 
the Balkans seems to have expanded to the Italian Peninsula. The 
land bridge which existed in the northern and central part of the 
present Adriatic Sea between the Italian and the Balkan 
peninsulas, during the Quaternary cold periods [6], would have 



made the contact between these two peninsulas easier. The 
expansion was followed by a slight differentiation in Italy from the 
Balkans. The COI haplotype network shows that the dispersal of 
haplotypes from the Iberian Peninsula to the north of Italy has also 
occurred, maybe either crossing mountains or along the Mediter- 
ranean coast of Spain and France as suggested for Cicada orni [52]. 
Although a weak flier, it is also possible that P. spumarius dispersed 
over the sea facilitated by wind (anemohydrochoric dispersal) since 
this mode of dispersal has already been observed in this species 
[20]. 

The expansion dates estimated here have wide confidence 
intervals. However, the lowest boundaries of these dates are about 
0.015 Ma ago, suggesting that the demographic and spatial 
expansion of this species may have occurred earlier than 
Holocene. Our data and reference [26] suggest the current 
geographic pattern of the species seems to result from postglacial 
range expansion from the Iberian Peninsula to the central and 
north-western Europe and, from the Anatolia/Caucasus (and 
eventually from western Asia) to east, north and central Europe, 
thus seemingly not following one of the common four paradigms 
[12,53]. Although a northern expansion from Balkans cannot be 
completely ruled out, the current data indicates that the 
Carpathians may have represented a geographic barrier to the 
northern expansion of Balkans populations. Further detailed 
sampling and genetic analysis of the Carpathian region would 
be important to test this hypothesis. 

Contact zones in Europe have been recorded for several 
European temperate species (reviewed in [54]). According to our 
results, the mtDNA lineages are geographically separated in most 
part of the range of P. spumarius but came into contact in some 
geographic regions (Fig. 1). This suggests the existence of recent 
admixture (secondary contact of diverged lineages) between 
mtDNA lineages in populations from these regions, also corrob- 
orated by [26]. The presence of the "eastern-Mediterranean" sub- 
lineage in some populations from the eastern part of the Iberian 
Peninsula (haplotype H56) and of the "western" sub-lineage in 
Balkans (haplotype HI 8) indicates that recent migrations between 
Mediterranean refugia may have occurred during the Quaternary 
period as reported in the olive fly Bactrocera oleae [55]. However, 
incomplete lineage sorting of an ancestral polymorphism cannot 
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Table 3. Tajima's D and [44] Fu's Fs test values and their statistical significance for Philaenus spumarius Cytochrome c oxidase I 
mtDNA groups. 





Neutrality Tests 




Tajima's Dtest 


Fu's Fs test 


Eastern Group 


0.15553 


-6.28375** 


Western-Mediterranean Group 


-1.77941* 


-23.61561*** 


Eastern-Mediterranean Group 


-1.75709* 


-26.22826*** 


*: indicates significant values at P<0.05; ** 


: indicates significant values at P<0.01 and ***: indicates 


significant values at P<0.001. 



doi:1 0.1 371 /joumal.pone.0098375.t003 



be ruled out as another possible explanation for the current 
mtDNA pattern of P. spumarius. Our data also indicate the 
existence of incomplete lineage sorting and/ or admixture in the 
nuclear gene EF- 1 a. Although nuclear eludes A and B correspond 
well to the "eastern-Mediterranean" and to the "western- 
Mediterranean" mtDNA sub-lineages, respectively, clade C is a 
mix of individuals from the three mitochondrial lineages. 
Heterozygous individuals whose alleles grouped within different 
clades were found, a fact not detected by [26], since they analysed 
homozygous individuals only. In the nuclear gene there was a lack 
of support for the monophyly of the "western" lineage ("western- 
Mediterranean and "eastern-Mediterranean" mitochondrial sub- 
lineages). Taken together, our results suggest that the current 
biogeographic pattern of P. spumarius may be the result of both 
secondary admixture and incomplete lineage sorting. 

Also quite interesting is the uncommon [12] high genetic 
diversity detected in P. spumarius populations from northern 
Europe (Scandinavia) indicating that the north of Europe was 
colonized by populations that may have survived in several extra- 
Mediterranean glacial refugia in addition to the "classical" 
Mediterranean ones [7]. 

Gene-flow between the Iberian Peninsula and Morocco 

The presence of the "western-Mediterranean" sub-lineage in 
Morocco suggests a close relationship between these P. spumarius 
populations and the Iberian Peninsula. This close relationship is 
also corroborated by the nuclear data. There is evidence that the 
Strait of Gibraltar has not been an effective barrier to the 
dispersal, having been, in fact, the route of dispersal for many 
species from Africa to Europe and vice-versa [56,57,58]. Lowered 
sea levels during glacial periods possibly facilitated exchange 
across the Strait of Gibraltar [59] . It is quite possible that, during 
such lower sea level periods, individuals from the Iberian 
Peninsula reached North Africa via anemohydrochoric dispersal. 
Contrarily to the thermophilous species Cicada barbara, a common 
cicada in southern Portugal and Spain [60] , the Rif Mountains did 
not appear to have acted as a geographical barrier to the dispersal 
of P. spumarius through Morocco, since this latter can be found in a 
variety of terrestrial habitats and even at altitudes above 1 700m 
(e.g. Mt Parnassus, in Greece: observations by Ana Rodrigues, 
Sara Silva and Eduardo Marabuto). Although the haplotypes 
found in samples from Morocco were the same as some of the ones 
found in the Iberian Peninsula, indicating that they belong to P. 
spumarius, the analysis of four male genitalia from these populations 
revealed similarities with P. tesselatus [61] and showed the necessity 
of further investigation on the taxonomic status of these species, as 
previously suggested by [26]. 



The UK and the origin of the North American and insular 
populations 

Our analyses suggest the presence of at least two mitochondrial 
lineages in the UK, the "western-Mediterranean" and the 
"eastern", and support a British origin for the populations of the 
Azores and New Zealand, and a multiple origin for the North 
American populations (Figs. 1 and 2). Populations from the Azores 
and New Zealand seem to have originated from only one of the 
British lineages here represented (the "western-Mediterranean" 
lineage). North American populations seem to have a mixed 
origin: a British origin, already suggested in the preliminary study 
by [27], from the "western-Mediterranean" lineage present in the 
UK, and an Iberian origin, due to the close relationship between 
one haplotype (H26) from the eastern United States and the 
Iberian Peninsula haplotypes. In fact, multiple translocations from 
different localities from western Europe have already been 
suggested for North American populations [26]. The close 
relationship of some North American haplotypes (New Hamp- 
shire) to the ones found in Anatolia and Finland indicates another 
origin from the "eastern" mitochondrial lineage, and that was 
never reported before. Nevertheless, a morphological variation in 
North American populations was already reported by [62]. The 
author shows that P. spumarius populations from New Hampshire 
and adjacent areas of North American exhibit morphological 
variation in male genitalia features and attributes this variation to 
hybridization between P. spumarius subspecies from different parts 
of Europe. Verification of whether there is any correlation 
between Hamilton's morphological subspecies and haplotype 
variation would require a parallel investigation of morphology 
beyond the scope of this work. The colonization of New Zealand 
and North America was probably recent and resulted from non- 
intentional anthropogenic introductions. This recent colonization 
could explain the spitdebug pest status in some locations of the 
USA and Canada, where P. spumarius reaches high densities, 
perhaps as result of the lack of competitors and predators [19,24]. 
Concerning the Azorean populations, the fact that they only occur 
in high altitude native habitats (e.g. in the geologically older areas 
of S. Miguel), very far from human altered habitats, is somewhat 
puzzling and we cannot exclude the possibility of a natural 
colonization by long-distance dispersal. 

Conclusion 

P. spumarius is one of the most widespread insects in Europe. We 
successfully provided time estimates of the main demographic and 
evolutionary events for the populations occurring in the main 
Mediterranean peninsulas, and in addition interpreted the 
colonization patterns out of Eurasia, namely of north-western 
Africa, North America, and the islands of Azores and New 
Zealand. This combination of a well analysed phylogeographic 
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and demographic pattern with the multiple transcontinental 
colonization events, some putatively natural, others recent non- 
intentional anthropogenic introductions, with ecosystem level 
consequences, make this species well placed for understanding 
the long term effects of invasive species and their post-invasion 
evolution. 

Supporting Information 

Figure SI Median-joining haplotype network of a set of 
Philaenus spumarius sampled geographic regions for 
mitochondrial gene COII (495bp). Size of the circles is in 
proportion to the number of haplotypes. Branches begin in the 
centre of the circles and their size is in proportion to the number of 
mutations. 
(TIF) 

Figure S2 Median-joining haplotype network of a set of 
P. spumarius sampled geographic regions for mitochon- 
drial gene cyt b (434bp) . Size of the circles is in proportion to 
the number of haplotypes. Branches begin in the centre of the 
circles and their size is in proportion to the number of mutations. 
(TIF) 

Figure S3 Median-joining haplotype network of P. 
spumarius sampled geographic regions for mitochon- 
drial gene COI (289bp). Size of the circles is in proportion to 
the number of haplotypes. Branches begin in the centre of the 
circles and their size is in proportion to the number of mutations. 
(TIF) 

Figure S4 Mismatch distribution of mtDNA COI haplo- 
types for each of the three P. spumarius haplogroups. (a) 

Eastern lineage; (b) Western-Mediterranean sub-lineage and (c) 
Eastern-Mediterranean sub-lineage. The expected frequency is 
based on a population growth-decline model, determined using 
DnaSP and is represented by a continuous line. The observed 
frequency is represented by a dotted line. Parameter values for the 
mismatch distribution are given in Table 2. 
(TIF) 

Figure S5 Bayesian skyline plots showing the historical 
demographic trends for each main Philaenus spumarius 
mtDNA group detected using COI gene. Along the y-axis is 
the expressed population size estimated in units of Neu, (Ne: 
effective population size, u.: mutation rate per haplotype per 
generation). The y-axis is in a log-scale. Solid lines represent 
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